rm(list=ls())
load("../Data/revealedLegis.RData")
set.seed(73964)
adapt <- 1e4
burn <- 25e4
samp <- 1e4
thin <- 25
chains <- 3

group.counter <- unique(allmeps$party)
group <- rep(NA,nrow(allmeps))
for (i in 1:length(group.counter)){
  group<- ifelse(allmeps$party==group.counter[i],i,group)
}
n.groups <- length(levels(as.factor(allmeps$party)))

nat.counter <- unique(allmeps$state)
nat <- rep(NA,nrow(allmeps))
for (i in 1:length(nat.counter)){
  nat<- ifelse(allmeps$state==nat.counter[i],i,nat)
}
n.nat <- length(levels(as.factor(allmeps$state)))

year <- allmeps$ep -3
n.year <- length(unique(year))

allmeps$nationalopen <- ifelse(is.na(allmeps$nationalopen),0,allmeps$nationalopen)
allmeps$nationalclosed <- ifelse(is.na(allmeps$nationalclosed),0,allmeps$nationalclosed)
allmeps$euopen <- ifelse(is.na(allmeps$euopen),0,allmeps$euopen)
allmeps$natback <- ifelse(is.na(allmeps$natback),0,allmeps$natback)
allmeps$other <- ifelse(is.na(allmeps$other),0,allmeps$other)
allmeps$epback <- ifelse(is.na(allmeps$epback),0,allmeps$epback)
allmeps$groupRole <- ifelse(is.na(allmeps$groupRole),0,allmeps$groupRole)
allmeps$committeeRole <- ifelse(is.na(allmeps$committeeRole),0,allmeps$committeeRole)

X <- with(allmeps,cbind(nationalopen,nationalclosed,euopen,epback,natback,other,age,groupRole,committeeRole))
n.betas <- ncol(X)
Jagsdata <- list(r=allmeps$votes,n=allmeps$totalvotes,n.obs=nrow(allmeps),
                 X=X,group=group,n.groups=n.groups,n.nat=n.nat,nat=nat,
                 year=year,n.year=n.year,n.betas=n.betas,
                 b0=rep(0,n.betas),B0=diag(.001,n.betas))


initsfunction <- function(chain) return(switch(chain,
                                               "1"=list(deltaTmp=rnorm(n.groups,0,2),gammaTmp=rnorm(n.nat,1,2),beta=rnorm(n.betas,0,3),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Super-Duper",.RNG.seed=2),
                                               "2"=list(deltaTmp=rnorm(n.groups,-1,3),gammaTmp=rnorm(n.nat,-2,3),beta=rnorm(n.betas,-1,2),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Wichmann-Hill",.RNG.seed=2),
                                               "3"=list(deltaTmp=rnorm(n.groups,1,2),gammaTmp=rnorm(n.nat,0,3),beta=rnorm(n.betas,1,1.5),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Mersenne-Twister",.RNG.seed=2)))

foo <- run.jags(model="../jags/binomialhierSurv.jags",
                data=Jagsdata,inits=initsfunction,modules="glm",
                adapt=adapt,burnin=burn,sample=samp,
                n.chains=chains,method='parallel',thin=thin,#summarise=FALSE,plots=FALSE,
                monitor=c("beta","gamma","deviance","delta","nu"))

save(foo,file="../Posteriors/ResultsVotesLegis.RData")

#### Days
allmeps2 <- subset(allmeps,subset=days>0)
X <- with(allmeps2,cbind(nationalopen,nationalclosed,euopen,epback,natback,other,age,groupRole,committeeRole))
n.betas <- ncol(X)
Jagsdata <- list(r=allmeps2$days,n=allmeps2$totaldays,n.obs=nrow(allmeps2),
                 X=X,group=group,n.groups=n.groups,n.nat=n.nat,nat=nat,
                 year=year,n.year=n.year,n.betas=n.betas,
                 b0=rep(0,n.betas),B0=diag(.001,n.betas))


initsfunction <- function(chain) return(switch(chain,
                                               "1"=list(deltaTmp=rnorm(n.groups,0,2),gammaTmp=rnorm(n.nat,1,2),beta=rnorm(n.betas,0,3),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Super-Duper",.RNG.seed=2), 
                                               "2"=list(deltaTmp=rnorm(n.groups,-1,3),gammaTmp=rnorm(n.nat,-2,3),beta=rnorm(n.betas,-1,2),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Wichmann-Hill",.RNG.seed=2),
                                               "3"=list(deltaTmp=rnorm(n.groups,1,2),gammaTmp=rnorm(n.nat,0,3),beta=rnorm(n.betas,1,1.5),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Mersenne-Twister",.RNG.seed=2)))

foo <- run.jags(model="../jags/binomialhier.jags",
                data=Jagsdata,inits=initsfunction,modules="glm",
                adapt=adapt,burnin=burn,sample=samp,
                n.chains=chains,method='parallel',thin=thin,#summarise=FALSE,plots=FALSE,
                monitor=c("beta","gamma","deviance","delta","nu"))

save(foo,file="../Posteriors/ResultsVotesLegisDays.RData")

#### Electoral systems
X <- with(allmeps,cbind(natCLPR,euCLPR,euCLPR_STV,natSemi_OLPR,
                        euSemi_OLPR,natSTV,euSTV,euOLPR,natSMP,epback,natback,
                        other,age,groupRole,committeeRole))
n.betas <- ncol(X)
Jagsdata <- list(r=allmeps$votes,n=allmeps$totalvotes,n.obs=nrow(allmeps),
                 X=X,group=group,n.groups=n.groups,n.nat=n.nat,nat=nat,
                 year=year,n.year=n.year,n.betas=n.betas,
                 b0=rep(0,n.betas),B0=diag(.001,n.betas))

initsfunction <- function(chain) return(switch(chain,
                                               "1"=list(deltaTmp=rnorm(n.groups,0,2),gammaTmp=rnorm(n.nat,1,2),beta=rnorm(n.betas,0,3),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Super-Duper",.RNG.seed=2),
                                               "2"=list(deltaTmp=rnorm(n.groups,-1,3),gammaTmp=rnorm(n.nat,-2,3),beta=rnorm(n.betas,-1,2),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Wichmann-Hill",.RNG.seed=2),
                                               "3"=list(deltaTmp=rnorm(n.groups,1,2),gammaTmp=rnorm(n.nat,0,3),beta=rnorm(n.betas,1,1.5),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Mersenne-Twister",.RNG.seed=2)))

foo <- run.jags(model="../jags/binomialhierSystem.jags",
                data=Jagsdata,inits=initsfunction,modules="glm",
                adapt=adapt,burnin=burn,sample=samp,
                n.chains=chains,method='parallel',thin=thin,#summarise=FALSE,plots=FALSE,
                monitor=c("beta","gamma","deviance","delta","nu"))

save(foo,file="../Posteriors/ResultsVotesLegisSystem.RData")

load("../Data/surveyLegis.RData")
set.seed(873209)

group.counter <- unique(outdata$party)
group <- rep(NA,nrow(outdata))
for (i in 1:length(group.counter)){
  group<- ifelse(outdata$party==group.counter[i],i,group)
}
n.groups <- length(levels(as.factor(outdata$party)))
colnames(outdata)[6] <- "state"
nat.counter <- unique(outdata$state)
nat <- rep(NA,nrow(outdata))
for (i in 1:length(nat.counter)){
  nat<- ifelse(outdata$state==nat.counter[i],i,nat)
}
n.nat <- length(levels(as.factor(outdata$state)))

year <- outdata$ep -4
n.year <- length(unique(year))

outdata$nationalopen <- ifelse(is.na(outdata$nationalopen),0,outdata$nationalopen)
outdata$nationalclosed <- ifelse(is.na(outdata$nationalclosed),0,outdata$nationalclosed)
outdata$euopen <- ifelse(is.na(outdata$euopen),0,outdata$euopen)
outdata$natback <- ifelse(is.na(outdata$natback),0,outdata$natback)
outdata$other <- ifelse(is.na(outdata$other),0,outdata$other)
outdata$epback <- ifelse(is.na(outdata$epback),0,outdata$epback)
outdata$groupRole <- ifelse(is.na(outdata$groupRole),0,outdata$groupRole)
outdata$committeeRole <- ifelse(is.na(outdata$committeeRole),0,outdata$committeeRole)

X <- with(outdata,cbind(nationalopen,nationalclosed,euopen,epback,natback,other,age,groupRole,committeeRole))
n.betas <- ncol(X)
Jagsdata <- list(r=outdata$votes,n=outdata$totalvotes,n.obs=nrow(outdata),
                 X=X,group=group,n.groups=n.groups,n.nat=n.nat,nat=nat,
                 year=year,n.year=n.year,n.betas=n.betas,
                 b0=rep(0,n.betas),B0=diag(.001,n.betas))


initsfunction <- function(chain) return(switch(chain,
                                               "1"=list(deltaTmp=rnorm(n.groups,0,2),gammaTmp=rnorm(n.nat,1,2),beta=rnorm(n.betas,0,3),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Super-Duper",.RNG.seed=2), 
                                               "2"=list(deltaTmp=rnorm(n.groups,-1,3),gammaTmp=rnorm(n.nat,-2,3),beta=rnorm(n.betas,-1,2),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Wichmann-Hill",.RNG.seed=2),
                                               "3"=list(deltaTmp=rnorm(n.groups,1,2),gammaTmp=rnorm(n.nat,0,3),beta=rnorm(n.betas,1,1.5),
                                                        nuTmp=rnorm(n.year),
                                                        .RNG.name="base::Mersenne-Twister",.RNG.seed=2)))

foo <- run.jags(model="../jags/binomialhier.jags",
                data=Jagsdata,inits=initsfunction,modules="glm",
                adapt=adapt,burnin=burn,sample=samp,
                n.chains=chains,method='parallel',thin=thin,#summarise=FALSE,plots=FALSE,
                monitor=c("beta","gamma","deviance","delta","nu"))

save(foo,file="../Posteriors/ResultsVotesLegisSurv.RData")